function [bo] = convert_rhombohedral_to_hexagonal(bo)
%CONVERT_RHOMBOHEDRAL_TO_HEXAGONAL Summary of this function goes here
%   Detailed explanation goes here

rx = bo.BasisVectors(:, 1);
ry = bo.BasisVectors(:, 2);
rz = bo.BasisVectors(:, 3);

hx = (2*rx-ry-rz)/3;
hy = (rx+ry-2*rz)/3;
hz = (rx+ry+rz)/3;

basisVectors = [hx hy hz];
lengthBasisVectors = size(basisVectors, 1);
add = [0 0 0;
  2/3 1/3 1/3;
1/3 2/3 2/3];

bo.BasisVectors = zeros(lengthBasisVectors*size(add, 1), 3);

for i = 1:size(add, 1)
  bo.BasisVectors(lengthBasisVectors*(i-1)+1:i*lengthBasisVectors, :) = ...
    basisVectors + repmat(add(i, :), lengthBasisVectors, 1);
end

ind = find(bo.BasisVectors > 1);
bo.BasisVectors(ind) = bo.BasisVectors(ind) - 1;

ind = find(bo.BasisVectors < 0);
bo.BasisVectors(ind) = bo.BasisVectors(ind) + 1;

bo.Atoms = [bo.Atoms; bo.Atoms; bo.Atoms];

[bo.BasisVectors, bo.Atoms] = sort_atoms_by_type(bo.BasisVectors, bo.Atoms);

end

